#ifndef _SPLINE_H_
#define _SPLINE_H_

#include <iostream>
#include <limits>
#include <cmath>
#include <algorithm>
#include <vector>
#include "function.h"
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/Sparse>
using namespace std;
using namespace Eigen;

class Interpolation
{
public:
    virtual void solve() = 0; 
    virtual double operator()(double _x){
        return 0;  
    }  
};

class ppFormInterpolation : public Interpolation
{
private:
    int n, condition, ord; 
    vector<double> X;
    Function &f;
    vector<Polynomial> P; 
    vector<vector<double>> C;
    VectorXd A_m;
public:
    ppFormInterpolation(Function &_f, vector<double> _X, int _condition, int _ord): f(_f), X(_X), condition(_condition), ord(_ord){
        n = X.size();
    }
    void solve()
    {
        //解线性方程组
        if (ord == 3)
        {
            SparseMatrix<double> A(n, n);
            vector<Triplet<double>> T;
            MatrixXd y(n, 1); 
            y = MatrixXd::Zero(n, 1);
            for (int i = 1; i < n - 1; i++)
		    {
                T.push_back(Triplet<double>(i, i - 1, (X[i + 1] - X[i]) / (X[i + 1] - X[i - 1])));
                T.push_back(Triplet<double>(i, i, 2.0));
                T.push_back(Triplet<double>(i, i + 1, (X[i] - X[i - 1]) / (X[i + 1] - X[i - 1])));
                y(i, 0) = 3.0 * (X[i] - X[i - 1]) / (X[i + 1] - X[i - 1]) * (f(X[i + 1]) - f(X[i])) / (X[i + 1] - X[i]) + 3 * (X[i + 1] - X[i]) / (X[i + 1] - X[i - 1]) * (f(X[i]) - f(X[i - 1])) / (X[i] - X[i - 1]);
            }

            if (condition == 1)
            {
                T.push_back(Triplet<double>(0, 0, 1.0));
                T.push_back(Triplet<double>(n - 1, n - 1, 1.0));
                y(0, 0) = f.diff1(X[0]);
                y(n - 1, 0) = f.diff1(X[n - 1]);
            }
            else
            {
                T.push_back(Triplet<double>(0, 0, 4.0));
                T.push_back(Triplet<double>(0, 1, 2.0));
                T.push_back(Triplet<double>(n - 1, n - 2, 2.0));
                T.push_back(Triplet<double>(n - 1, n - 1, 4.0));
                if (condition ==  2)
                {
                    y(0, 0) = 6.0 * (f(X[1]) - f(X[0])) / (X[1] - X[0]) - f.diff2(X[0]) * (X[1] - X[0]);
                    y(n - 1, 0) = 6.0 * (f(X[n - 1]) - f(X[n - 2])) / (X[n - 1] - X[n - 2]) + f.diff2(X[n - 1]) * (X[n - 1] - X[n - 2]);
                }
                if (condition ==  3)
                {
                    y(0, 0) = 6.0 * (f(X[1]) - f(X[0])) / (X[1] - X[0]);
                    y(n - 1, 0) = 6.0 * (f(X[n - 1]) - f(X[n - 2])) / (X[n - 1] - X[n - 2]);
                }
            }

            A.setFromTriplets(T.begin(), T.end());
            A.makeCompressed();
            SparseLU<SparseMatrix<double>> S;
            S.compute(A);
            A_m = S.solve(y);

            for (int i = 1; i < n; i++)
		    {
                vector<double> pp;
                double K;
                K = (f(X[i]) - f(X[i - 1])) / (X[i] - X[i - 1]);
                pp.push_back(f(X[i - 1]));
                pp.push_back(A_m(i - 1));
                pp.push_back((3.0 * K - 2.0 * A_m(i - 1) - A_m(i)) / (X[i] - X[i - 1]));
                pp.push_back((A_m(i - 1) + A_m(i) - 2.0 * K) / pow(X[i] - X[i - 1], 2));
                C.push_back(pp);
                Polynomial p(pp, X[i - 1]);
                P.push_back(p);
            }
        }
        if (ord == 1)
        {
            for (int i = 1; i < n; i++)
		    {
                vector<double> pp;
                double K;
                K = (f(X[i]) - f(X[i - 1])) / (X[i] - X[i - 1]);
                pp.push_back(f(X[i]) - K * X[i]);
                pp.push_back(K);
                C.push_back(pp);
                Polynomial p(pp, 0);
                P.push_back(p); 
		    }
        }
    }
    double operator()(double _x)
    {
        if ((_x < X[0]) || (_x > X[n - 1]))
        {
            cout<< "Error!" <<endl;
		    exit(0);
        }
        if (_x == X[0])
            return f(X[0]);
        int point = 1;
        while (_x > X[point])
        {
            point = point + 1;
        }

        if (_x == X[point])
            return f(X[point]);
        else
            return P[point-1](_x);
    }
    vector<vector<double>> Get_A()
    {
        return C;
    }
};

class BsplineInterpolation : public Interpolation
{
private:
    int n, condition, ord; 
    vector<double> X, x1, x2;
    vector<Bsplines> B;
    Function &f;
    VectorXd A_x; 
public:
    BsplineInterpolation(Function &_f, vector<double> _X, int _condition, int _ord): f(_f), X(_X), condition(_condition), ord(_ord)
    {
        n = X.size();
        if (ord == 1)
        {
            x1 = {X[0] - 1 * (X[1] - X[0])};
            x2 = {X[n - 1] + 1 * (X[1] - X[0])};
        }
        if (ord == 3)
        {
            x1 = {X[0] - 3 * (X[1] - X[0]), X[0] - 2 * (X[1] - X[0]), X[0] - 1 * (X[1] - X[0])};
            x2 = {X[n - 1] + 1 * (X[1] - X[0]), X[n - 1] + 2 * (X[1] - X[0]), X[n - 1] + 3 * (X[1] - X[0])};
        }
        X.insert(X.begin(), x1.begin(), x1.end());
        X.insert(X.end(), x2.begin(), x2.end());
        for (int i = 2 - ord; i < n + 1; i++){
            vector<double> v(X.begin() + i - 2 + ord, X.begin() + i + 2 * ord);
            Bsplines B1(ord, v);
            B.push_back(B1);
        }
    }
    double operator()(double _x)
    {
        if (ord == 1)
        {
            if (_x == X[1])
                return f(X[1]);

            int point = 1;
            while (_x > X[point + 1])
            {
                point = point + 1;
            }
            return A_x(point - 1) * B[point - 1](_x) + A_x(point) * B[point](_x);
        }
        if (ord == 3)
        {
            if (_x == X[3])
                return f(X[3]);
            int point = 1;
            while (_x > X[point+3])
            {
                point = point + 1;
            }
            return A_x(point - 1) * B[point - 1](_x) + A_x(point) * B[point](_x) + A_x(point + 1) * B[point + 1](_x) + A_x(point + 2) * B[point + 2](_x);
        }
        return 0;    
    }
    void solve()
    {
        //求解线性方程组
        if (ord == 1)
        {
            SparseMatrix<double> A(n, n);
            vector<Triplet<double>> T;
            MatrixXd y(n, 1); 
            y = MatrixXd::Zero(n, 1);
            for (int i = 0; i < n; i++){
                T.push_back(Triplet<double>(i, i, 1.0));
                y(i, 0) = f(X[i + 1]);
            }
            A.setFromTriplets(T.begin(), T.end());
            A.makeCompressed();
            SparseLU<SparseMatrix<double>> S;
            S.compute(A);
            A_x = S.solve(y);
        }
        if (ord == 3)
        {
            SparseMatrix<double> A(n+2, n+2);
            vector<Triplet<double>> T;
            MatrixXd y(n+2,1); 
            y = MatrixXd::Zero(n+2, 1);
            for (int i = 1; i < n + 1; i++){
                T.push_back(Triplet<double>(i, i - 1, B[i - 1](X[i + 2])));
                T.push_back(Triplet<double>(i, i, B[i](X[i + 2])));
                T.push_back(Triplet<double>(i, i + 1, B[i + 1](X[i + 2])));
                y(i, 0) = f(X[i + 2]);
            }

            if (condition ==  1)
            {
                vector<double> v0(X.begin() + 1, X.begin() + 5);
                vector<double> v1(X.begin() + 2, X.begin() + 6);
                vector<double> vn_1(X.end() - 6, X.end() - 2);
                vector<double> vn(X.end() - 5, X.end() - 1);
                Bsplines B0(2, v0), B1(2, v1), Bn1(2, vn_1), Bn(2, vn);
                T.push_back(Triplet<double>(0, 0, (-3.0 * B0(X[3])) / (X[4] - X[1])));
                T.push_back(Triplet<double>(0, 1, (3.0 * B0(X[3])) / (X[4] - X[1]) - (3.0 * B1(X[3])) / (X[5] - X[2])));
                T.push_back(Triplet<double>(0, 2, (3.0 * B1(X[3])) / (X[5] - X[2])));
                T.push_back(Triplet<double>(n + 1, n - 1, (-3.0 * Bn1(X[n + 2])) / (X[n + 3] - X[n])));
                T.push_back(Triplet<double>(n + 1, n, (3.0 * Bn1(X[n + 2])) / (X[n + 3] - X[n]) - (3.0 * Bn(X[n + 2])) / (X[n + 4] - X[n + 1])));
                T.push_back(Triplet<double>(n + 1, n + 1, (3.0 * Bn(X[n + 2])) / (X[n + 4] - X[n + 1])));
                y(0, 0) = f.diff1(X[3]);
                y(n+1, 0) = f.diff1(X[n+2]);
            }
            else
            {
                T.push_back(Triplet<double>(0, 0, 6.0 / ((X[4] - X[1]) * (X[4] - X[2]))));
                T.push_back(Triplet<double>(0, 1, -6.0 / ((X[4] - X[1]) * (X[4] - X[2])) - 6.0 / ((X[5] - X[2]) * (X[4] - X[2]))));
                T.push_back(Triplet<double>(0, 2, 6.0 / ((X[5] - X[2]) * (X[4] - X[2]))));
                T.push_back(Triplet<double>(n + 1, n - 1, 6.0 / ((X[n + 3] - X[n]) * (X[n + 3] - X[n + 1]))));
                T.push_back(Triplet<double>(n + 1, n, -6.0 / ((X[n + 3] - X[n]) * (X[n + 3] - X[n + 1])) - 6.0 / ((X[n + 4] - X[n + 1]) * (X[n + 3] - X[n + 1]))));
                T.push_back(Triplet<double>(n + 1, n + 1, 6.0 / ((X[n + 4] - X[n + 1]) * (X[n + 3] - X[n + 1]))));
                if (condition ==  2)
                {
                    y(0, 0) = f.diff2(X[3]);
                    y(n + 1, 0) = f.diff2(X[n + 2]);
                }
                if (condition == 3)
                {
                    y(0, 0) = 0.0;
                    y(n + 1, 0) = 0.0;
                }
            }

            A.setFromTriplets(T.begin(), T.end());
            A.makeCompressed();
            SparseLU<SparseMatrix<double> > S;
            S.compute(A);
            A_x = S.solve(y);
        }
    }
    VectorXd Get_A()
    {
        return A_x;
    }
};

class Curvespline : public Interpolation
{
private:
    int n, d, m, ord; 
    const double interval = 0.01;
    vector<vector<double> > X, Fit_X;
    vector<double> ccl;
public:
    Curvespline(vector<vector<double> > _X, int _m, int _ord): X(_X), m(_m), ord(_ord)
    {
        n = X.size();
        d = X[0].size();
        for (int i = 0; i < n; i++){
            if (X[i].size() != d){
                cout << "Error!" << endl;
                exit(0);
            }
		}
    }
    double norm2(vector<double> x, vector<double> y)
    {
        double dist = 0.0;
        for (int i = 0; i < d; i++){
            dist = dist + (x[i] - y[i]) * (x[i] - y[i]);
        }
        return sqrt(dist);
    }
    void solve(){
        ccl.push_back(0.0);
        for (int i = 1; i < n; i++)
		{
            ccl.push_back(ccl[i - 1] + norm2(X[i], X[i - 1]));
        }
        vector<double> E(n);
        for (int i = 0; i < d; i++)
		{
            vector<double> temp1, temp2;
            for (int j = 0; j < n; j++)
		    {
                temp1.push_back(X[j][i]);
            }
            Discrete temp_f(ccl, temp1, E, E);
            if (m == 1)
            {
                BsplineInterpolation BSI(temp_f, ccl, 3, ord);
                BSI.solve();
                double l = 0;
                while (l <= ccl[n - 1])
                {
                    temp2.push_back(BSI(l));
                    l = l + interval;
                }
            }
            if (m == 2)
            {
                ppFormInterpolation ppF(temp_f, ccl, 3, ord);
                ppF.solve();
                double l = 0;
                while (l <= ccl[n - 1])
                {
                    temp2.push_back(ppF(l));
                    l = l + interval;
                }
            }
            Fit_X.push_back(temp2);
        }
    }
    vector<vector<double> > Get_X()
    {
        return Fit_X;
    }
};


#endif
